A posteriori error approximation in discontinuous Galerkin method on polygonal meshes in elliptic problems

The paper presents a posteriori error approximation concept based on residuals in the two-dimensional discontinuous Galerkin (DG) method. The approach is relatively simple and effective in application, and it takes advantage of some unique properties of the DG method. The error function is constructed in an enriched approximation space, utilizing the hierarchical nature of the basis functions. Among many versions of the DG method, the most popular one is based on the interior penalty approach. However, in this paper a DG method with finite difference (DGFD) is utilized, where the continuity of the approximate solution is enforced by finite difference conditions applied on the mesh skeleton. In the DG methods arbitrarily shaped finite elements can be used, so in this paper the meshes with polygonal finite elements are considered, including quadrilateral and triangular elements. Some benchmark examples are presented, in which Poisson’s and linear elasticity problems are considered. The examples use various mesh densities and approximation orders to evaluate the errors. The error estimation maps, generated for the discussed tests, indicate a good correlation with the exact errors. In the last example, the error approximation concept is applied for an adaptive hp mesh refinement.

This paper presents an error approximation concept for a version of the discontinuous Galerkin (DG) method called the discontinuous Galerkin with finite difference (DGFD) method, introduced in the first author's previous papers [1][2][3][4][5] . The main characteristic of the DGFD method is a stabilization parameter w , interpreted as a small distance from the mesh skeleton. In the DGFD method, like in any other version of the DG method, arbitrary global approximation can be constructed on the entire domain with arbitrary basis functions in elements. In this paper Chebyshev polynomials are employed as the approximation basis since their construction does not require nodes and, as such, they do not depend on the shape of finite elements. In consequence, arbitrary polygonal finite elements can be used in the mesh. The Chebyshev polynomials can be constructed in a recursive manner, which helps one to avoid significant truncation errors when constructing higher-order polynomials. These basis function properties are utilized for error approximation.
Error estimation is a common problem in all computational methods. It is used to assess the quality of the approximate solution. When it is applied for mesh refinement, error estimation makes it possible to obtain a better solution with minimum computational costs. In any automatic mesh adaptation, the procedure of error estimation is the most important component. Therefore, this paper is focused on error estimation techniques in the DGFD method to be applied in an automatic hp mesh refinement.
The DG method has a long history that goes back to the early seventies of the previous century 6 . A constant development has been observed since that time [7][8][9] . In general, the DG method can be treated as a particular version of the finite element method. In the DG method the domain is discretized using a finite element mesh, and the basis functions are constructed within the finite elements. The main difference between those two methods is the construction of the global approximation space. In the standard FEM, the global approximation space is continuous a-priori, while in the DG method the global approximation space is not continuous, so some techniques have to be employed to make the final approximate solution continuous. The lack of such continuity in the DG method seems to be its drawback, but on the other hand, it provides great flexibility in applying the basis functions and the shapes of the finite elements. It is natural for the DG methods to use polygonal/polyhedral finite elements. Moreover, the orders of approximation in the finite elements can be selected independently for each

Problem formulation
This paper utilizes two types of two-dimensional linear elliptic boundary value problems to present a new error approximation method. The problems are defined on domain Ω with outer boundary Γ . The first one is the scalar Poisson's problem with the following boundary conditions where Γ D and Γ N are the parts of the outer boundary with Dirichlet and Neumann boundary conditions, respectively, and Γ D ∪ Γ N = Γ , Γ D ∩ Γ N = ∅.
The second problem is the two-dimensional plane strain elasticity problem which can be expressed in the following way where σ is the stress tensor, b is the body force vector, t is the prescribed traction force vector, u is the displacement vector and û is the prescribed displacement vector.
Equation (2) is supplemented by Hooke's law where ε is the strain tensor, E is the fourth order Hooke's tensor, in which for the two-dimensional plane strain state the following non-zero terms are present and , µ are standard elasticity constants. Small strains are assumed in this paper, therefore the Cauchy strain tensor is used In the DGFD method the discontinuities of the global approximation across the mesh skeleton Γ d need to be considered in the problem formulation. Local coordinates associated with the mesh skeleton are specified. They are based on two unit vectors (n, s) , where n is the unit vector normal to Γ d and s is tangent to Γ d . Fig. 1a presents the skeleton local coordinates, and Fig. 1b shows a graphical illustration of a stability parameter w . The parameter w is related only to the mesh, not to the considered problem. Therefore, different types of boundary problems can be solved using the same mesh with a similar value of the stability parameter 63 . The numerical fluxes on the mesh skeleton are evaluated using finite difference formulas in the small neighborhood of the mesh skeleton which is defined by parameter w 3 . In the DGFD method, w is not only a penalty-like parameter. It indicates points in the vicinity of the mesh skeleton, see Fig. 1b, where the trial function or its derivatives should be calculated to evaluate the numerical fluxes. It has been shown in other papers, for example [3][4][5]63 , that the DGFD method is stable and not sensitive to fluctuations of w . A wider discussion of parameter w is presented in "Discussion on the stability parameter" Section. In the discontinuous Galerkin method, the discontinuity and mean-value operators appear, which are defined with the help of the vector n normal to the mesh skeleton. Their definitions for an auxiliary function g are as follows:

Weak formulations
Following the procedures given in 2,5 the problems presented in Eqs. (1) and (2) can be rewritten in weak forms. In both cases of the weak forms the test and trial functions belong to the same broken Sobolev space S , defined as follows where Ω h represents the set of disjoint finite elements covering the whole domain Ω , and Ω e is a single finite element. In the space S a subset S ⊂ S can be distinguished that consists of functions which are continuous on the mesh skeleton and meet the boundary values on the Dirichlet boundary In this section the weak form of Poisson's problem is first presented, and then, avoiding certain repetitions, the weak form of the elasticity problem is shown. The weak form of Poisson's problem in Eq. (1) is expressed as follows: Find u ∈S such that where the bilinear and linear forms are defined as It is assumed in Eq. (9) that �∇u� · n = 0 on Γ d , which is a standard assumption in the computational methods such as FEM or other versions of the DG method. In this weak form, the test function belongs to the space of discontinuous functions, which is inconvenient when we want to find the approximate solution, since we want the trial function to belong to S . Therefore, this formulation must be modified so that both the test and trial functions belong to the same space, however, the trial functions are forced to belong to set S when the stability parameter w tends to zero. The modified weak form of Poisson's problem is now expressed as follows:Find u ∈ S such that: where where x w = x − w n , x 2w = x − 2 w n and the value of the trial function at point x w on the Dirichlet boundary means the value of this function at distance w from the boundary in the normal direction: u(x w ) = u(x − nw) for x ∈ Γ D . Now, we present the lemma showing that the formulations in Eqs. (9) and (11) are equivalent.

Lemma 1
If there is a function u ∈ S that satisfies Eq. (11), then the function belongs to S and satisfies Eq. (9).
Proof Equation (11) has to be satisfied for every value of w , and, in particular for an arbitrarily small value of the parameter. This means that every finite difference relation in this equation has to tend to a specific derivative for w going to zero. Therefore, on the mesh skeleton we have ∀ v ∈ S and for w → 0 www.nature.com/scientificreports/ The finite difference on the mesh skeleton has to be finite for an arbitrarily small value of w , which indicates that the trial function has to be continuous on the mesh skeleton. On the Dirichlet boundary, we have the following finite difference relation In this case the finite difference on the Dirchlet boundary is also finite for every value of w , which indicates that the trial function has to meet the boundary values on the Dirichlet boundary. Additionally, it is obvious that the mean values of the derivative on the mesh skeleton tend to the derivative on the mesh skeleton, and a similar situation occurs with the derivative on the Dirichlet boundary Equation (11) is now rewritten with some small modifications When w goes to zero, then the finite difference relations in the second and fourth integrals in Eq. (16) tend to the derivatives as shown in Eqs. (13) and (14), and, using relations (15), this results in the following equation in which u ∈S . After a few easy algebraic operations in Eq. (17), the Eq. (9) is retrieved.
To obtain a numerical solution, Eq. (11) has to be rewritten in a discrete version, in which the infinite functional space S is substituted by the finite space S p ⊂ S . The subspace S p comprises a complete set of polynomial basis functions in each element up to order p . In computations parameter w cannot be arbitrarily small, the parameter is just small in the numerical sense. The Poisson's problem shown in Eq. (11) is expressed in the discrete form as follows:Find u p ∈ S p in Ω such that The boundary value problem of elasticity in Eq. (2) also has its discrete form obtained in the same way as for the Poisson's problem. It was derived in 4,5 and reads Find u p ∈ S p in Ω such that The elasticity problem is a vector problem, so the definitions of the bilinear and linear forms are more complex: The tensors E d i and E b i , related to the mesh skeleton and outer boundary, respectively, are defined as follows ∀ v p ∈ S p and for small value of w

Approximation in DG method
In the DG methods, the basis functions in the elements can be chosen as any sort of polynomial basis functions.
In particular, they can just be monomials, Chebyshev, Legendre, or any other sort of polynomials, provided that they are complete up to order p . In certain situations, the basis functions can be enriched by another sort of functions, see 3 . In this paper, the approximation basis is chosen to be the Chebyshev polynomials, which are constructed using the recursive form as follows 64 Orthogonality is a well-known property of the Chebyshev polynomials. However, they are not purely orthogonal, but orthogonal with the weight 1 and so on. The global approximation, for the entire domain Ω , is constructed with the help of the global approximation matrix p and the vector of degrees of freedom L u . The suffix p in the symbol of the approximation matrix means that the polynomials of order up to p are used in the elements to compute the matrix. For example, the approximations for the scalar field u and the vector field u are as follows The matrix p is a one-row matrix for the scalar function u in the Poisson's problem given in Eq. (1) and two-row matrix for the 2D elasticity problem in Eq. (2).
In the standard finite element method, the approximation is constructed with the help of the so-called shape functions, whose definitions depend on the shape of the finite element. In the DG methods, the basis functions do not depend on the shapes of the elements. The element matrices in the DG methods have the hierarchic structure as illustrated in Fig. 2. This means that the stiffness matrix of order p keeps inside all the stiffness matrices of the lower orders. On the other hand, when we have the stiffness matrix of order p , it is relatively easy to extend it to the matrix of order p + 1.
When the approximation in Eq. (26) is applied to the weak formulation of the problem in Eq. (11) or (18), we obtain the discrete problem in an algebraic form where K p and f p are defined depending on the problem under consideration  www.nature.com/scientificreports/ The DGFD method has been programmed in the Matlab environment which offers great flexibility and many up-to-date mathematical tools. Particularly, the subroutines for linear solvers are available in this environment. In this paper, the 'linsolve' function has been applied for solving algebraic equations. The main problem in Eq. (27) is solved in the standard way.

Discussion on the stability parameter
In Eqs. (18) and (19) the scalar stabilisation parameter w plays an important role. It should cause the final discrete solution to be continuous and the Dirichlet boundary conditions to be met. The value of the parameter should be small enough for the compatibility conditions to be satisfied, but not too small to avoid numerical instabilities.
The w parameter measures the small distance from the mesh skeleton, so its value depends on the size of the adjacent elements. On the other hand, its value is independent of the considered problem, but has to be adjusted to the mesh, so it is possible to use a single stabilization parameter for coupled problems 63 . The parameter cannot be arbitrarily small since numerical instabilities can occur if it is too small. In consequence, the solutions obtained this way are not strictly continuous and the Dirichlet boundary conditions are not exactly met. In other words, the value of the parameter should be chosen small enough to obtain good quality solutions, but not too small to avoid truncation errors that can spoil the solution.
The w parameter should be small, however, it turns out that the final solution is not sensitive to its value. This means that a relatively large change in the value causes tiny changes in the approximate solution. This indicates that the DGFD method is consistent and well-defined. In the DGFD method, elements of various sizes and orders can be freely combined in the mesh, so the method is especially eligible for mesh adaptation. The stabilization parameter depends only on the mesh, i.e. its value depends only on the size of elements that adhere to the skeleton segment, w = γ · h , where h indicates the element size and γ is a scaling parameter that is constant for the whole mesh skeleton.
In order to illustrate the influence of the stabilization parameter on the final solution in the DGFD method, the Poisson's problem, see Eq. (1), is solved in the square domain Ω = [−1, 1] × [−1, 1] . The right-hand side function and the boundary conditions are taken from the exact solution, which is assumed to be the following exponential function with the parameters α = 10 and β = 2 . The function with its both derivatives is depicted in Fig. 3.
The problem has been solved with various values of scaling parameter γ for quadrilateral and polygonal meshes. The meshes and maps of errors are shown in Figs. 4 and 5, respectively. In the quadrilateral mesh, the largest element is over one thousand times larger than the smallest element, and the two elements are adjacent to each other. The polynomial orders in the two meshes have been randomly generated: p ∈ [3,20] for the quadrilateral mesh and p ∈ [3,8] for the polygonal mesh. It can be seen in the two examples that even though the value of the stabilization parameter is changed by five orders of magnitude, only slight differences in the errors can be observed. Thus, the stability parameter in the DGFD method can be used to combine finite elements of different sizes as well as different orders.

Error estimation method
In the proposed method of error approximation, the original boundary value problem is considered in a residuum form. The weak form of an error function is constructed, which is used to obtain an approximate error function. For the sake of clarity, the method is first presented for Poisson's problem and then extended to the elasticity problem. The derivation begins from defining an exact error function as the difference between the exact solution u and its approximation obtained by solving Eq. (18) in Ω (32) e =ê =û − u p on Γ D , ∇e · n =ĥ − ∇u p · n on Γ N  The test as well as the trial functions in Eq. (33) belong to the infinite space, that is v, e ∈ S , while the approximate solution belongs to the discrete space u p ∈ S p . We want to find an approximation of error e , and so the problem defined in Eq. (33) is rewritten in the finite space as follows:For given function u p ∈ S p find approximation ep ∈ Sp such that The symbol p indicates the order of approximation space, in which the error function is approximated. It is possible that p is the same as p , but much better results are obtained when p > p . In this paper, it is assumed that p = p + 1 . As it is noticed in "Approximation in DG method" Section, the enriched space in the DGFD method is quite easy to construct due to the hierarchical form of the basis functions. A similar error approximation equation to the one presented in Eq. (34) can be generated for the elasticity problem. In that case the final problem of error approximation reads:For given u p find ep ∈ Sp in Ω such that where The error function ep is approximated in the same manner as in Eq. (26), and using Eqs. (34) or (35a) the following system of algebraic equations is obtained ∀ v ∈ S and for w → 0 ∀ vp ∈ Sp and for small value of w  www.nature.com/scientificreports/ The approximation in the elements is constructed in the hierarchic way, as shown in "Approximation in DG method" Section and so, and therefore the system of algebraic equations in (36) is arranged in the following way The matrix K p has already been used to solve the Eq. (27), so it is already factorized. This is utilized in the error approximation The solution scheme shown in Eq. (38) requires the solution of the relatively small algebraic problem to get vector L e 2 and then vector L e 1 is obtained. The reduced system of algebraic Eq. (38a) is connected with the enriched part of the approximation space, which is much smaller in comparison to the main problem and therefore it is especially effective for high-order approximation.

Examples
This section presents four benchmark examples, illustrating the performance of the developed error estimation method. Poisson's problem is considered in the first two examples ("Exponential example for Poisson's problem, Hyperbolic example for Poisson's problem" Sections). The third example focuses on an elasticity problem ("Elasticity problem with singularity" Section). In the fourth example the Poisson's problem is considered once again and an hp mesh adaptation is performed using the error estimation approach to identify the elements for the refinement. Error estimates are compared in each case with the exact error for various mesh densities and approximation orders. In order to present the quality of the error estimation the following efficiency index is defined for the approximate error in the whole domain and in a single finite element Exponential example for Poisson's problem. In this section, a benchmark Poisson's problem is discussed whose exact solution is known a priori. It is the exponential function shown in Eq. (29) with α = 10 and β = 2 , depicted in Fig. 6a. It can be observed that the function is non-zero in the region close to the domain center, thus, the error concentrations are expected to appear there. Three polygonal meshes shown in Fig. 6b- Tables 1, 2 and 3 where the global error measure e L 2 (Ω) is compared with the approximate global error �ep� L 2 (Ω) with p = p + 1 and the efficiency indices are also provided in these tables. It is shown that the value of �ep� L 2 (Ω) gets closer and closer to the exact value with the increase of mesh density and approximation order. It can be observed in Tables 1, 2 and 3 that the best results have been obtained for the polygonal meshes, while the least accurate results are for the triangular mesh. The rate of convergence for all kinds of meshes is the same, however, the error estimates for the triangular meshes show some irregularity, especially for higher order approximation. In the triangular elements, the vertex inner angles are acute, which may lead to some perturbations when the continuity of the approximate solution is enforced. It is the reason why the error estimation for the triangular elements does not show a stable performance, especially for high-order approximation. The in Ω → KpL e = rp www.nature.com/scientificreports/ polygonal finite elements, i.e. quadrilateral, pentagonal, and hexagonal elements, usually give better results in comparison to the triangular ones, which has been confirmed in this example. For a deeper analysis, the results have been presented in the form of error maps for three meshes: mp250, msq256, and mtr242 as well as for three values of p = 3, 5, 8 . The results are shown in Figs. 7, 8 and 9 in the form of maps of exact and approximate errors for the three meshes and the three approximation orders. As it is seen in the figures, the exact error is quite well recovered by the approximate error ep . In each case, the error distributions are reconstructed with a relatively high accuracy using the proposed method.
Tab. 1 presents the global efficiency indexes for the polygonal meshes. For a deeper analysis of the error, the maps of the efficiency indexes over the fifth order elements for three meshes m50, m25 and m1000 are shown in Fig. 10. Comparing the maps of efficiency indexes with the error maps in Fig. 7 it can be concluded that the efficiency indexes are very close to one for the elements located at the places of error concentration. At the places where the error level is very small the element efficiency index tilts from unity. This means that the error approximation technique works well at the places with error concentration, however, at the places where the level of the error is very small its approximation may differ from the exact error values.
So far, the error approximation methods discussed in this paper have been applied only to meshes with uniform approximation orders. However, the method is also suitable for unstructured meshes with non-uniform    www.nature.com/scientificreports/ elements orders. The Poisson's exponential benchmark has also been solved using a quadrilateral, randomly refined mesh where the orders of the elements have been chosen randomly from the range p = 2 to p = 10 . The mesh structure and the map of the orders of the elements are depicted in Fig. 11. The analyzed error approximation technique is able to cope with such a strongly heterogeneous mesh because the error can be approximated with a relatively good accuracy as shown in Fig. 12.     www.nature.com/scientificreports/ The domain has been discretized with the polygonal mesh mp1000, see Fig. 6d. Uniform approximation orders p = 3, 5, 8 have been applied to all elements in the mesh. The global errors, their approximations and efficiency indices are shown in Table 4. The approximate global errors, once again, are very close to the exact errors. The maps of errors and their approximations are presented in Fig. 14, showing accurate concentrations of errors produced by the new error approximation technique.
The solution of Eq. (41) is singular at the origin since the stresses at this point tend to infinity. That is why in the approximate solution the error concentration is expected to appear around this point. The benchmark was originally solved on the so-called rotated L-shaped domain, see 65 . In this paper, the standard L-shaped domain with vertices at points: (− 1, − 1), (1, − 1), (1, 0), (0, 0), (0, 1), (− 1, 1) is considered. To obtain the solution in such a domain, both the global coordinates and the displacements need to be appropriately transformed, as shown below. The polar coordinates are constructed using Cartesian coordinates ξ and η where the auxiliary coordinates (ξ , η) come from the following transformation of the global coordinates  The maps of errors for the polygonal mesh, provided in Fig. 16a,b, enable the comparison of the exact and approximated errors. In this example, the error is concentrated, as expected, in the vicinity of the origin. The error approximation procedure is able to locate the place with the error concentrations and correctly recover the level of the error. The global efficiency index in this example is ηp = 0.5 , which can be caused by the error concentration in the very narrow region in the vicinity of the origin. For a deeper analysis of the situation the maps of the element efficiency indexes are presented in Fig. 16c. It can be noticed that in a large part of the domain the element efficiency indexes are close to one in most of the elements except those located in two regions where the indexes are about 0.7, however, the regions are not located near the origin. This example confirms that the procedure for error approximation performs well also for the elasticity problem.  In order to obtain an accurate solution by means of the discrete method the mesh should be refined along the narrow band. In this example, an automatic hp mesh refinement is applied using the a posteriori error approximation presented in this paper. In this example quadrilateral elements are used, since they enable flexible refinement in the DGFD method, as shown in the example in "Exponential example for Poisson's problem" Section. In each step of the refinement, the problem is solved on the current mesh with the error approximation. Then, the error level for each element is calculated.
Various scenarios for the combined h and p mesh refinement can be applied. The main role of the automatic mesh adaptation is to minimize the global error and make the error uniformly distributed in the domain. The mesh adaptation is performed in refinement steps, in which the error is estimated after the problem solution on the current mesh. It is important to prepare a proper refinement scenario to reach the final mesh with a relatively small number of refinement steps. In this application the error is estimated for every finite element, then the elements with maximum and minimum errors are identified. In every adaptation step, two error parameters are set: the first is the 0.2 and 0.8 combination of the maximum and minimum errors in elements, respectively, and the second is the 0.05 and 0.95 combination of the same errors. The two error parameters are applied in the hp mesh refinement. The finite elements in which the error is higher than the second error parameter are p refined by increasing the approximation order by two and they are divided into four elements. Moreover, their neighbors are p refined by one and h refined by dividing them into four elements. The unrefined elements with errors greater than the first error parameter are divided into four elements. The whole procedure starts with a very coarse mesh, 3 × 3 with second-order finite elements. In the first three steps, only the h refinement is applied to locate the places with higher errors. Afterward, the combined hp refinement is applied as described above. Many hp schemes can be applied to reach a similar final refined mesh. However, in this case, the p and h refinement are applied in every step of the refinement procedure to limit the number of iterations. Using this approach we obtained the final refined mesh in relatively few steps. It should be emphasized that we got the refined mesh that consists of elements with a wide spectrum of orders, from p = 2 up to p = 13 . In the refinement procedure, the large and small finite elements as well as the elements with low and high order of approximation can be neighbors. It shows that the properties of the DGFD method allow for very flexible mesh refinement. The refinement of the p or h type can be applied only to a single element in the DGFD method without the need to interfere with the neighboring elements. In this example, the ability of the DGFD method to perform effectively combined mesh refinement is presented, and automated adaptivity is based on the effective error estimation in the elements.
Selected meshes obtained during the refinement procedure are shown in Fig. 18. It is seen that the mesh is mainly refined along the band where the errors are concentrated, so that the global error can be reduced with a relatively small number of degrees of freedom ( #dofs ). In the first stage of the process, only the h-refinement is performed in order to identify the places with error concentration. In Fig. 19a the convergence diagram for the refinement is depicted. This is visible in the quite flat first part of the convergence diagram. After this stage, further refinement is practically concentrated in the vicinity of the specific band, which results in the steep second part of the diagram. The value of the efficiency index ηp during the refinement procedure is shown in Fig. 19b. For the first five steps of the refinement the efficiency index is underestimated, however in the further steps of the refinement procedure the index values are close to one. The initial perturbations in the efficiency index are caused by the coarse mesh structure and the fact that the error is concentrated in the narrow band in the domain. For the refined mesh the value of the global error is correctly estimated. In the fifth step the efficiency index is about -8, where the real error was 1.5 × 10 0 and the approximated error was 4.2 × 10 −2 . In Fig. 20 the maps of the element efficiency indexes are presented for the fifth, ninth and thirteenth adaptation steps. It can be noticed   www.nature.com/scientificreports/ that in each case the indexes are close to one near the band with error concentration and the element indexes with higher values than one are generally at places with small error values.

Conclusions
The paper presents the method for error approximation in the DGFD method for elliptic problems analyzed with various polygonal meshes. Some unique properties of the DGFD method are utilized to obtain an effective error approximation in the approximate solution. As a result, a new technique for error approximation is proposed. Chebyshev polynomials are applied as the basis functions in the DGFD method. Their definition is recursive, thus they are hierarchical, which enables an efficient generation of the matrices for higher-order approximation space. The hierarchical formation of the stiffness matrix is utilized to get the error approximation developed in the upgraded approximation space.
The stability parameter w plays an important role in the DGFD method, since it is used for evaluating the socalled numerical fluxes on the mesh skeleton and for enforcing the continuity of the final solution. Generally, the stability parameter should be small for the compatibility of the final solution. However, the DGFD method is not very sensitive to the variations of this parameter, i.e. it can be changed significantly with a very small influence on the final solution. This property is utilized in this paper for error approximation techniques.
The method of error approximation has been illustrated with four benchmark examples, in which two kinds of two-dimensional elliptic problems, Poisson's and elasticity problem, have been analyzed. It has been shown that the error approximation method on the polygonal meshes is able to locate places with error concentrations. In the last example, the method has been applied with success for an automatic hp mesh refinement. The method of error approximation can recover the exact error measure of high quality for polygonal meshes, including quadrilateral and triangular ones.
So far the error approximation method for the DGFD method has been applied in 2D elliptic problems. In the further research the method will be applied to other problems like linear and non-linear Navier-Stokes problem, or diffusion-convection transport in a compressible or incompressible medium. www.nature.com/scientificreports/